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Abstract. - The recent experimental realization of cold polar molecules in the rotational and 
vibrational ground state opens the door to the study of a wealth of phenomena involving long- 
range interactions. By applying an optical lattice to a gas of cold polar molecules one can create 
a layered system of planar traps. Due to the long-range dipole-dipole interaction one expects a 
rich structure of bound complexes in this geometry. We study the bilayer case and determine the 
two-body bound state properties as a function of the interaction strength. The results clearly 
show that a least one bound state will always be present in the system. In addition, bound states 
at zero energy show universal behavior and extend to very large radii. These results suggest 
that non-trivial bound complexes of more than two particles are likely in the bilayer and in more 
complicated chain structures in multi-layer systems. 



> 

oo 
in 

cn 
o 
o 



X 



Introduction. — Quantum gases of polar atoms and 
molecules in their rovibrational ground-state represent a 
unique opportunity to study the interplay of long- and 
short-range interactions in the highly controllable trapped 
gas environment. Early experiments used magnetic dipo- 
lar atoms [IH1] which have observable effects in spite of 
intrinsically weak dipole moments. Recently, heteronu- 
clear molecules with very large electric dipole moments 
have been realized by a number of groups [SHPJ ■ The goal 
of a quantum degenerate system of polar molecules with 
strong 1/r 3 long-range dipole-dipole forces therefore seems 
close at hand. 

The attractive force of polar molecules in the head- 
to-tail configuration can lead to collapse of the system 
[ID] . However, as suggested by Wang et al. [IT], a one- 
dimensional optical lattice that creates a multilayered 
stack of pancake systems can stabilize the situation. If 
we apply a field to polarize the dipoles perpendicular to 
the layers then the intralayer interaction will be purely re- 
pulsive, whereas the interlayer part will be attractive but 
with the optical lattices separating the dipoles in different 
layers. This setup is presently being implemented exper- 
imentally. As discussed in [IT], the dipole-dipole force 
forms bound chains and the system effectively behaves as 
a liquid of chains with resemblance to rheological fluids. 



In the case of bosons we expect the chains to Bose con- 
dense at low temperatures. However, if we have fermionic 
polar molecules the situation is less clear since one would 
expect a Bose-Fermi mixtures with various bound com- 
plexes [12] . The dipole potential can also now be inverted 
to have a repulsive core with the use of laser fields [T5irH] , 
which offers different physics possibilities than the "natu- 
ral" dipole orientation. 

If we simplify the problem to consider just two adjacent 
layers we have a system with 1 /r 3 interactions that mim- 
ics the long-range 1/r interactions in graphene |15| and 
semiconductor bilayers [16) . In the semiconductor case 
bound states of exciton pairs with non-zero dipole mo- 
ments have been considered in connection with organic 
interfaces and quantum wells [17] . For small coupling 
strength it was concluded that no bound state exists \T7\ . 
This was also stated in several recent works concerning 
cold polar molecules [UJ[T21[TH] where the conclusion was 
based on a Gaussian ansatz. However, at small coupling 
the particles are strongly delocalized and a localized Gaus- 
sian is therefore not a good approximation [19]. To make 
matters worse, the potential integrates to zero over the 
plane and thus at small coupling the Landau criterion [20] 
for a bound state in two-dimensional systems is not appli- 
cable. Using scattering theory it was recently shown that 
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Fig. 1: Potentials of different shapes, dipole-dipole potential 
(Vdip), t wo square wells (V a tt, V rep ), harmonic oscillator (Vho), 
and f/r (Vi/ r ). All potentials have the nodes at r = R s , the 
same attractive and repulsive volumes, and net volume zero. 

a bound state presumably exists for arbitrarily small mo- 
ments 21,22 . However, the scattering theory is intricate 
and does not yield straightforward information about the 
behavior of the two-body bound state wave function. 

The purpose of the present work is to compute and ex- 
plain the basic properties of two-body systems used as 
the fundamental building blocks for layered dipolar struc- 
tures. We shall employ simple model potentials to extract 
universal properties, point out where details of the po- 
tentials are needed, and illustrate relations between wave 
functions, energies, and radii. We shall use square well, 
harmonic oscillator and 1/r-potentials in two dimensions, 
and compare to solutions of the true dipole-dipole poten- 
tial. We briefly sketch the model solutions, discuss ener- 
gies, threshold properties, and various implications. 

Model solutions. — We solve the 2-dimensional (2D) 
Schrodinger equation, which is possible to do analyti- 
cally with several model potentials. We use cylindrical 
coordinates (r, 9) and separate the total wave function 
* = R(r)G>(0) into radial R(r) and angular <J>(#) parts. 
With Q(9) = exp(im6')/v / 27r, where m is an integer, the 
stationary radial Schrodinger equation becomes: 

-2 



2m 



dr 2 



1 - 4m' 
4r 2 



V(r) 



u(r) = Eu(r), (1) 



where u — ^frR(r) is the reduced radial wave function, 
E is the energy, and V(r) is the potential, assumed to 
be spherically symmetric. We consider potentials with 
net volume equal to zero, J V(r)d 2 r = 0, just as for the 
dipole-dipole potential, i.e. 



V(r) = D 2 



-2d 2 



(r 2 + d 2 fl 2 



(2) 



where D is the electric dipole moment and d is the distance 
between the two different layers containing the particles. 



This potential is in Fig. [TJ compared with square well, har- 
monic oscillator, and 1/r potentials at distances less than 
their common node at R s = d^/2. The equal volume con- 
ditions relate strengths and shifts for given radial shapes. 
Thus we have D 2 = V S R% (3/2) 3 / 2 /2, and for the square 
well shape, the small, V s , and large, V\, distance absolute 
strength values are related by 



M5r> ' 



(3) 



where Ri is the radius where the outer square well ends. 
The attractive part is the most interesting and we shall 
only use the harmonic and 1 Jr potentials with the square 
barrier shape for Ri > r > R s as illustrated in Fig. [TJ 

One of our main concerns is the appearance of a bound 
state at small couplings. We therefore only consider the 
most attractive m = potential which simplifies to 



dr 2 



A+ fc2 



u = 0, 



(4) 



where k 2 — 2m(E — V(r))/% 2 depends on r. 

Let us now first solve completely the piecewise con- 
stant potential in Fig. Q] for use as a reference standard. 
Then the wave functions in the three different regions of 
space are the Riccati Bessel functions of order —1/2, which 
means that the solutions for R(r) are various Bessel func- 
tions of order depending on the region of space: 



AJo(kir) r<R s 
R = { BI (k 2 r) + CK (k 2 r) R s < r < Ri 
DK (k 3 r) r>Ri, 



(5) 



where Jq,Iq,Kq are ordinary and modified Bessel func- 
tions of the first and second kind, respectively. The coeffi- 
cients A, B, C, and D are determined by matching at the 
region boundaries and by normalization. The wavenum- 
bers, fcj, are the absolute values of k in the regions. Match- 
ing logarithmic derivatives leads to the transcendental 
equation for the energies: 



k2RiK a (k 3 Ri)K 1 (k 2 Ri) - k 3 R l K (k 2 R l )K 1 (k 3 R l ) = 
k 3 RiI a (k 2 Ri)Ki(k 3 Ri) + k 2 RiK (k 3 Ri)Ii(k 2 Ri) 
fc 2 J? s J (fcii? s )^i(fc 2 i? s ) - fcifi s K (fc 2 i? s )Ji(fci J R s ) 
k 2 R s J Q (k 1 R s )I 1 (k 2 R s ) + k 1 R s Io(k 2 R s )J 1 (k 1 R s ) ' 



(6) 



This formula and the wave function, Eq.([5]), are valid 
for the potential V aU in Fig. [TJ For V rep from Fig. [JJ 
one merely takes the analytic continuation of the rele- 
vant Bessel function. Expanding the solution for small 
strengths results in 



E = E exp 



^£ (o + E ° 
'Vi \ VtHRt/Rs) 



(7) 



where Eo — 2h/(mRf). Eq.© with the second order po- 
tential strengths replaces the Landau expression when the 
potential integrates to zero |20) . 
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Fig. 2: Energies of states as a function of V s in units of 
h 2 /(mRs). Right and left correspond to V a tt and Vrep, re- 
spectively, of Fig. [T] The dotted grey lines are the calculated 
asymptotic values for the square well states. The solid red lines 
are the square well energies with Ri — 2R S , the dashed green 
lines are square well energies with Ri — 3R S , and the black 
double-dotted line the dipole-dipole potential results. The in- 
set is a zoom of small energy plotted logarithmically (values 
in points colored as in the main plot), the lines are from Eq. 
(J7J. The magenta and blue points, from 1/r and oscillator 
potentials, appear together with the dipole points found for 
-Ri = 3R 3 (Vi is from Eq. ((3}), and the green square well points. 
The black line for the dipole potential is to guide the eye. 



Analytic solutions can also be found when the short- 
distance part of the potential is substituted by a harmonic 
oscillator potential, I4 (r) = 2V s [(r/R s ) 2 — 1] and a 1/r 
potential, V\/ r {r) = V s (l—R s /r). In both cases the poten- 
tials are zero for r = R s and V s is the strength of the box 
potential with the same volume below R s . The radial wave 
functions for r < R s are given by R(r) — iVexp(— z/2)M, 
where M is the confluent hypergeometric function, or 
Kummer function M (a, b, z) with three arguments. We 
have b — 1 and a and z are: 



E + 2V S 
RsVs 



iRl 



Vs 



2r 2 y/mV 

hR s 



(8) 



in 



2{V S - E) ' 



for oscillator and 1/r potentials, respectively. This is the 
solution vanishing at r = 0. Matching at r — R s with the 
square well solutions for r > R s leads to trancendental 
equations for the energies. 

Energies. — The energies obtained from Eq. ([BJ are 
shown in Fig. [2] as function of V S R 2 , or equivalently of vol- 
ume of either attraction or repulsion. We choose Ri = 2R S 
or Ri = 3i2 S! as indicated, and use V from Eq.Q. In- 
creasing V s for both natural and inverted potentials leads 
to more bound states which more and more are determined 
from the attraction alone and independent of the confin- 



ing barriers. The condition, Jo(kiR s ) = 0, determines k\ 
for V a tt in Fig. Q]and thereby the bound state energy E n 
is related to the nodes of the Bessel function, i.e. 



E n 



-Vs 



ft a (*>, 



2mR 2 



(10) 



where j'q n is the nth zero of Jq(x). This asymptotic limit 
of straight lines of slope equal to one is valid for an in- 
finitely high barrier. The ground state is the deepest and 
agrees well with the asymptotic limit, but for the excited 
states the well must be deeper to reach the limit, around 
-150 for 1% agreement for the first excited state. 

For V rep in Fig. [TJ where the repulsion is for r < R s 
and the attractive well is at R s < r < Ri, no confining 
barrier exists at larger distance. Here, the bound state 
wave function is R(r) = AJq(t) + BNo(r), where No(r) is 
the Neumann function. In the deep well limit, the wave 
function must vanish at both endpoints where large energy 
implies large arguments of both the Bessel and Neumann 



functions. Then k^Ri — R S Y 



or equivalently 



E,, 



n 2 ir 2 h 2 



{Ri/R s ) 2 - 1 2m(R l - R s ) 2 ; 



(11) 



where the slopes of the lines are dependent on Ri and 
therefore different from the V a tt case. These estimates are 
much further off than the similar estimates for V a tt- For 
the difference between the energies and the estimate to be 
1% for the ground state, the depth of the well needs to 
be around 100. For that level of agreement in the first 
excited state, the well depth needs to be around 250. 

The result for small strengths are shown in the inset 
of FigO One bound state is always present for all po- 
tentials even when the strength is approaching zero. The 
dependence on both strengths and radius (Ri ) is substan- 
tially stronger for the inverted potentials. The very small 
energies close to threshold obey the limiting linear de- 
pendence from Eq. ([7]). They are extremely small due 
to the zero net volume, and hence difficult to obtain ac- 
curately. The numerical results for the dipole potential 
almost coincide with the energies for the other poten- 
tials. On the scale of the inset the energies are exceed- 
ingly small and difficult to calculate numerically, especially 
for the inverted dipole potential. For larger strengths 
mV s R 2 /h 2 > 2 on the inverted side and mV s R 2 /h 2 > 10 
on the natural side, we found that the energies rather pre- 
cisely are given by E = —coD 2 /d 3 exp(—cih 2 d/(mD 2 ), 
where (c ,ci) = (0.85, 9.86), (0.012, 19.9) for the natural 
and inverted potentials, respectively. 

The values in the experiments [B] correspond to 



strengths of about V s R 2 m/h — 1.9 or mR 2 E/h 



-0.6 



(Fig. [2]) which is on the verge of universality (see Fig. SJ, 
although current temperatures are too high to maintain 
such a weakly bound state [6] . Fortunately, systems with 
larger D 2 /d 3 can be explored where the bound state en- 
ergy is much larger and obeys the analytic formula given 
above. 
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Fig. 3: The lowest three bound state energies as function of 
strength for different potentials, square well (solid red), har- 
monic oscillator (dashed blue), 1/r (dotted magenta), and 
dipole (double-dotted black). The square well barrier is the 
same in all cases, Ri = 2R„. 



Fig. 4: Binding energy multiplied with < r 2 > plotted against 
the binding energy for the lowest bound states for the square 
well (parameters in the legend). The crosses are results for the 
dipole potential. The scale on the abscissa is logarithmic, with 
small binding energies to the left. 



These computations are also carried out for harmonic 
oscillator and 1/r potentials at small distance with a 
square well at larger distances. The condition of expo- 
nentially vanishing wavefunction with r implies that the 
first argument a of the Kummer function must be a non- 
positive integer — n counting the bound states. Then M 
reduces to the Laguerre polynomials. The energies become 



E. 



ho 



2[-V„ + 



E\/ r — V s - 



, rnRl 
2mV*R 2 s 
7L 2 (2n+l) 2 



(2n+l) 



(12) 
(13) 



where n can take the values 0,1,2, etc. These energies 
are only approached asymptotically as for the square well 
in Fig. [TJ They are highly potential dependent for well- 
bound states as the barrier essentially has no influence on 
these energies, see Fig.O The potential in Eq.fl2) leads to 
energies between those of harmonic and square well poten- 
tials. Thus limits to realistic potentials can be provided 
by analytic models. 

However, in the limit where V s — > both energies in 
Eqs. (TT2"]) and (IT3l) approach zero from positive values. The 
approximations are too crude, and we replace the r > 
R s potentials by the corresponding square well. We then 
find numerically that the energies for both these potentials 
always remain negative corresponding to bound states for 
all values of V s . Thus all the investigated potentials with 
zero net volume always have at least one bound state. 

Threshold properties. — The weakly bound states 
often reveal unique physics as for example Efimov and halo 
states with universal properties. The number of nodes 
as function of volume can be found rather precisely from 
Eqs. (|10H13"l) by solving for n with E„ = 0. The explicit 



result for harmonic oscillator, 2n + 1 = ^Jm\V s \Rl/h, 
whereas for the 1/r potential the right-hand side is a 
factor y/2 larger. Thus these potentials would have the 
same number of bound states if the volume of the oscilla- 
tor were twice as large as that of the 1/r-potential. The 
dipole-dipole potential seems from Fig. Q] to have interme- 
diate properties. However, the dependencies on strength 
for each energy differs substantially, see Eqs. (fT2l and fp~3|) . 

The same extrapolation to the threshold can be applied 
to the wave functions. We then find that the mean square 
radius is given as < r 2 /R? s >= 1/2 for the oscillator and 
5/8 to leading order for the 1/r potential. These thresh- 
old radii, obtained from asymptotic strong binding, are 
however qualitatively completely wrong when the energy 
is sufficiently close to zero. In Fig. [2] the straight lines ap- 
proaching zero bend over as the system attempts to stay 
bound for a weaker potential. The wave function is cor- 
respondingly leaking out under the barrier as the energy 
approaches zero. This is the effect producing nuclear ha- 
los |19j where the mean square radius for two particles 
in three dimensions becomes inversely propertional to the 
energy. The analogue here is that most of the probability 
is found for r > Ri where the wave function is Kq. This 
means that the mean square radius then approaches |23j 



< r z >= 



f~ a»\K (x)\ 2 dx 



i h z 



2m\E\ /°° x\K (x)\ 2 dx 3m\E\ 



(14) 



which is a universal result independent of both the par- 
ticular (s-wave) state considered and the shape of the at- 
tractive potential. 

The rate of approach to the asymptotic value in Eq. (fl~4| 
is seen in Figure 2] Large binding energies correspond 
to wave functions located in the attractive well. As 
the threshold is approached the wave function begins to 
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Fig. 5: The probability, J 3 |u| 2 dr, of being in the attractive 
part of the V a tt potential for the lowest three bound states 
shown as a function of (r 2 ). For V re p, all states behave as the 
ground states of both potentials. 

leak out and eventually ends up in the universal limit in 
Eq. (fTl| . However, how that limit is approached depends 
on the presence of a barrier. Without a barrier (V rep from 
Fig. [1]), all states, both ground and excited, approach the 
weak binding limit in the same way. This is precisely the 
halo effect [13] . 

The similar approach to universality is found for the 
ground state for V a tt from Fig. Q] since then the bar- 
rier vanishes with the binding energy. For excited states 
of Vatt, the barrier remains finite for vanishing binding 
energy. High excitation goes together with high barrier 
which causes even the weakly bound states to remain lo- 
calized in the attractive region. The approach to univer- 
sality is delayed by many orders of magnitude in binding 
energy. During the approach the radius remains small for 
small energies and consequently the curves dive below the 
universal limit before the eventual approach. These fea- 
tures become more pronounced with excitation energy. 

In Figure@]we also show results for different values of i?; , 
still with Vi adjusted to maintain zero net volume. The 
curves for the ground states are roughly identical. This 
is also seen at higher energies for the first excited state, 
but eventually at small energies a wider barrier delays the 
escape of the wave function even though the barrier is 
correspondingly smaller. 

The rate of approach to the universal limit at threshold 
is further illustrated in Figure [3] The probability distri- 
bution from the wavefunction is moving slower than the 
distribution in the mean square radius integral. The bulk 
of the probability can stay in the attractive well while the 
root mean square radius is much larger than the radius of 
the attraction. Thus the tail properties are crucial. These 
effects are enhanced for excited states because the final ex- 
pulsion of the wave function to the external region comes 
much later for states constraint by a large barrier. 



The formal connection is that the wave function essen- 
tially is Jo until the energy is very small where it has to 
change to Kq , see Eq. ([5]) . Computation of mean square ra- 
dius then employs Jo f° r normalization but Kq for the re- 
distribution. Thus neither short nor large-distance proper- 
ties are sufficient for a description in this transition region. 

Implications for many particles. Several conclu- 
sions are immediately deduced from our model systems. 
First of all, at least one bound state must appear in 2D 
for any potential with zero net volume. Thus, in contrast 
to the statements in jTT][T21[T71[IH], one bound state is 
always present in a bilayer with dipoles oriented perpen- 
dicular to the layers. This in turn means that arguments 
based on the existence of a critical strength for binding 
should be re-considered. Furthermore, since a large dipole 
moment implies a large strength, several two-body bound 
states will be present. This must be taken into account 
in simulations of configurations of actual systems where 
a finite temperature might lead to population of differ- 
ent excited states. In practice, finite temperature in the 
system will put a natural limit on how small binding en- 
ergies one needs to consider, i.e. for \E\ < ksT the bound 
states are thermally dissociated into the continuum and 
therefore largely irrelevant. 

Second, a positive net volume still allows bound states 
in two dimensions. However, now it is necessary to have 
a non-zero minimum attraction where the strength would 
increase with the net volume. This resembles the situation 
in three-dimensional quantum mechanics. 

Third, the structure of many particles in layers are 
strongly influenced by their pair interactions. The absence 
of a critical interaction strength for binding implies that a 
phase transition from a superfluid to a dipolar chain liquid 
cannot occur in contrast to the suggestion in [TTiri8] . Fur- 
thermore, the presence of bound states in bilayers for all 
strengths immediately implies that chains in multilayers 
will also be present for all couplings. We would therefore 
expect the system to always be a dipolar chain liquid. As 
all the individual chains can form immediately the chain- 
chain interaction becomes interesting and necesary to in- 
clude in careful future investigations. 

We can in fact give an upper bound on the binding 
energy of a chain of M particles in M different layers. The 
scaling of two-body energy with layer distance is E(nd) = 
Eijv?, with n an integer giving the distance in equally 
spaced layers a distance d apart. Here E2 is the two-body 
energy we have calculated above which depends on the 
coupling strength. Taking the large M limit, the upper 
bound is 

00 2 
E M =E 2 mY] -o = ^rM E 2 . (15) 

n— 1 

This estimate of the energy agrees well with the calcu- 
lation of |llj where the harmonic approximation was at 
small couplings for M = 2 and M = 81. 
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Fourth, the inverted potential with V rep in Fig. [T]with 
repulsion at small distances and an intermediate attractive 
pocket also has at least one bound state for any strength. 
In the bilayer system with a finite density of fermionic 
polar molecules in each layer, Ref. [T3] found an interesting 
particle-hole coherent state. The density used there was 
chosen small so intralayer interaction is negligible and the 
focus was on the interlayer repulsion corresponding to a 
barrier at small radii and an attractive pocket outside. 
In |13j the presence of bound states shown in this paper 
was ignored. In a real system at finite temperature this 
coherent state could probably exist at low dipole strength. 
However, at larger strengths the bound states must be 
taken into account and one would presumably have instead 
an interacting gas of bound pairs behaving as a superfluid. 

Fifth, the bound states become extremely extended as 
the energy gets sufficiently close to zero. Such delocalized 
two-body states can in turn enter into complicated multi- 
particle complexes [24]. A bound three-body system of 
two particles in one plane and one in between in the adja- 
cent plane is therefore possible. This presupposes that the 
repulsive intraplane interaction between these extended 
structures is less than two times the two-body binding en- 
ergy that created the two-body bound states. This sort of 
Y-junction configuration can be very interesting in ther- 
modynamic considerations of chains as it will contribute 
non-trivially to the entropy and can help lower the free 
energy. Full quantum studies of such configurations are 
therefore very relevant and worth pursuing. 

Sixth, the optimum conditions for interesting multi- 
particle structures are probably in the regime close to the 
threshold for binding of the second state. This can be be- 
low the threshold where additional attraction from other 
particles would lead to binding in analogy to Borromean 
three-body systems where the two-body subsystems are 
unbound [T9l[23] . It can also be for slightly larger attrac- 
tion and with a bound two-body state since such a system 
still is spatially extended and in a sense rather similar to 
the unbound continuum state. The latter case is analogous 
to the extremely weakly bound atomic helium dimer where 
the trimer becomes well bound but with a very weakly 
bound and spatially extended excited state [19ll23j . 

Seventh, the regime of weak binding, strong derealiza- 
tion and large root-mean-square radius exhibits universal 
features independent of the shape of the potentials. The 
same type of universality is likely to exist for multiple 
bound states but much finer tuning is probably required 
to reach these structures. 

Conclusions. — Presently experimentalists are work- 
ing to produce layered systems of dipolar molecules. We 
use simple model potentials to study the bilayer case with 
dipoles polarized perpendicular to the layers. We find the 
solutions and calculate properties of the wave functions 
and in particular binding energies and radii as functions 
of dipole moment or strength of the potentials. Realistic 
potentials are used to test the generality of our results. 



We conclude that there always is a bound state for all 
strengths of the dipolc-dipole potential. We find that the 
wave functions of both ground and excited states show 
universal behavior at zero energy as they basically reside 
where the potential has become vanishingly small. To 
access this universal regime the dipole strength must be 
tuned around the threshold for a bound state to appear. 
The extended wave functions indicate that three or more 
particle complexes are possible in chains and in bilayers. 
The repulsive in-plane interaction becomes interesting in 
connection with these structures. In any case the tuning 
of interactions to universal regimes emphasizes the close 
analogy to the physics studied through the well-known 
technique of Feshbach resonances. 



NTZ is grateful to B. Wunsch, D.-W. Wang and E. Dern- 
ier for numerous discussions. 
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